rm(list=ls())
# Set working directory 
setwd("")

# Load programs
source("procedures.R")
library(ggplot2)

# Random seed
set.seed(222)
# Sample size
n <- 1e+6
# Parameters to be saved here
thetas <- rep(NA, 6)
# Choose DGP
dgp.nr <- 4

# Graphics formatting parameters
mccolors <- function(...) gray.colors(..., end = 0.7) 
sel <- c(1, 3, 5, 6)
cols <- rep(c("gray70", "gray40", "gray20"), each = 2)[sel]
ltys <- rep(c(1, 2, 3), 2)[sel]
lwds <- c(3, 3, 6, 3, 3, 3)[sel]
nms <- c("Log", "Brier", "Spherical", "Boosting", "As1", "As2")[sel]
colb <- c(1, cols)
ltyb <- c(1, ltys)
lwdb <- c(5, lwds)
nmb <- c("True", nms)

# Simulate data
dat <- sim(dgp = dgp.nr, n = n)

# Load objective functions
of <- getof()

# Estimate parameter under various scoring rules
for (j in 1:6) thetas[j] <- optimize(of[[j]], y = dat$y, x = as.matrix(dat$x), interval = c(-10, 10))$minimum

# Draw prediction plot
pdat <- cbind(dat$p, sapply(thetas, function(s) logit(dat$grid*s)))
mccolors <- colorRampPalette(c("blue","red","green"))
pdf(paste0("pred", dgp.nr, "_", n, "_color.pdf"), family = "Bookman")
par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
matplot(x = dat$grid, y = pdat, type = "l", col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, 
        xlab = "x", ylab = "")
legend("topright", inset = c(-0.4, 0), c("True", "Log", "Brier", "Spherical", "Boosting", "As1", "As2"),
       col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, xpd = TRUE)
dev.off()

# Draw classification plot
pdf.options(family = "Bookman")
pdf(paste0("class", dgp.nr, "_", n, "_color.pdf"), family = "Bookman")
x <- 2
c <- as.matrix(seq(from=0.01, to=0.99, by=0.01))
tr <- (dat$p.fct(x) > c)
gd <- sapply(logit(x*thetas), function(s) s > c)
par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
matplot(cbind(tr, gd),type = "l", lwd=4, x=c, col=c(1,mccolors(6)), ylab="", xlab=expression(c), 
        lty=c(rep(1,4),rep(2,3)))
legend("topright", inset = c(-0.4, 0), c("True", "Log", "Brier", "Spherical", "Boosting", "As1", "As2"),
       col = c(1, mccolors(6)), lty = c(rep(1,4), rep(2,3)), lwd = 4, xpd = TRUE)
dev.off()

# Draw plot for share of correct predictions 
cgrid <- seq(from = 0.01, to = 0.99, by = 0.01)
x <- dat$x
true <- sapply(x, function(z) dat$p.fct(z) > cgrid)
share.correct <- matrix(0, length(cgrid), 6)
for (jj in 1:6){
  aux <- sapply(x, function(z) logit(z*thetas[jj]) > cgrid)
  share.correct[, jj] <- rowMeans(true == aux)
}

pdf(paste0("sc", dgp.nr, "_", n, "_color.pdf"), family = "Bookman")
par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
matplot(share.correct, type = "l", lwd = 4, x = cgrid, col=c(mccolors(6)), ylab="", 
        xlab=expression(c), lty=c(rep(1,3),rep(2,3)))
legend("topright", inset = c(-0.4, 0), c("Log", "Brier", "Spherical", "Boosting", "As1", "As2"),
       col = c(mccolors(6)), lty = c(rep(1,3), rep(2,3)), lwd = 4, xpd = TRUE)
dev.off()     